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We investigate a novel approach to measuring the Hubble constant using gravitational- wave (GW) 
signals from compact binaries by exploiting the narrowness of the distribution of masses of the 
underlying neutron-star population. Gravitational-wave observations with a network of detectors 
will permit a direct, independent measurement of the distance to the source systems. If the redshift 
of the source is known, these inspiraling double-neutron-star binary systems can be used as standard 
sirens to extract cosmological information. Unfortunately, the redshift and the system chirp mass 
are degenerate in GW observations. Thus, most previous work has assumed that the source redshift 
is obtained from electromagnetic counterparts. However, we investigate a novel method of using 
these systems as standard sirens with GW observations alone. In this paper, we explore what we 
can learn about the background cosmology and the mass distribution of neutron stars from the 
set of neutron-star (NS) mergers detected by such a network. We use a Bayesian formalism to 
analyze catalogs of NS-NS inspiral detections. We find that it is possible to constrain the Hubble 
constant, Ho, and the parameters of the NS mass function using gravitational-wave data alone, 
without relying on electromagnetic counterparts. Under reasonable assumptions, we will be able to 
determine Ho to ±10% using ~100 observations, provided the Gaussian half- width of the underlying 
double NS mass distribution is less than 0.04Mq. The expected precision depends linearly on the 
intrinsic width of the NS mass function, but has only a weak dependence on Ho near the default 
parameter values. Finally, we consider what happens if, for some fraction of our data catalog, we have 
an electromagnetically measured redshift. The detection, and cataloging, of these compact-object 
mergers will allow precision astronomy, and provide a determination of Ho which is independent of 
the local distance scale. 

PACS numbers: 98.80.Es, 04.30.Tv, 04.80.Nn, 95.85.Sz 
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I. INTRODUCTION 



The previous decade has seen several ground-based 
gravitational-wave (GW) interferometers built, and 
brought to their design sensitivity. The construction 
of Initial LIGO, the Laser Interferometer Gravitational- 
wave Observatory [l|, , was a key step in the quest for a 
direct detection of gravitational waves, which are a fun- 
damental prediction of Einstein's theory of gravity 0, 0] ■ 
The three LIGO detectors are located in the USA, with 
two sited in Hanford, Washington within a common vac- 
uum envelope (HI, H2 of arm- lengths 4 km and 2 km re- 
spectively) and one in Livingston, Louisiana (LI of arm- 
length 4 km) [HQ. The 600 m arm-length GEO-600 de- 
tector [5[ is located near Hannover, Germany. LIGO and 
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GEO-600 began science runs in 2002, and LIGO reached 
its initial design sensitivity in 2005. The 3 km Virgo 
interferometer Q, located at Cascina, Italy, began com- 
missioning runs in 2005, and has participated in joint 
searches with LIGO and GEO-600 since 2007. The 300 
m arm-length TAMA-300 detector [7|, located in Tokyo, 
Japan had undertaken nine observation runs by 2004 
to develop technologies for the proposed underground, 
cryogenically-cooled, 3 km arm- length LCGT project [8|. 

Gravitational waves from the coalescences of compact- 
object binaries [9| consisting of neutron stars (NSs) and 
black holes (BHs) are among the most promising sources 
for LIGO [l0( . The first joint search for compact binary 
coalescence signals using the LIGO S5 science run and 
the Virgo VSR1 data has not resulted in direct detec- 
tions, and the upper limits placed on the local NS-NS 
merger rate are higher than existing astrophysical upper 
limits 0- However, construction has already begun on 
the Advanced LIGO detectors llj, which are expected 
to increase the horizon distance for NS-NS inspirals from 
~33 to ~445 Mpc. This thousandfold increase in de- 
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tection volume is expected to yield detections of NS-NS 
coalescences at a rate between once per few years and 
several hundred per year, with a likely value of ^40 de- 
tections per year 12j. 

The Advanced Virgo detector [l3[ is expected to be- 
come operational on a similar timescale as Advanced 
LIGO (~2015) and with similar sensitivity. We denote 
the network of three AdLIGO detectors and AdVirgo as 
HHLV in the following. These may later be joined by 
additional detectors, such as LIGO Australia, IndlGO 
or LCGT, creating a world-wide detector network whose 
sensitivity will enable gravitational- wave astronomy. The 
network comprising LIGO-Australia, HI and LI will be 
denoted as AHL. Regardless of the precise network con- 
figuration, the hope is that when AdLIGO (and cither 
LIGO Australia or AdVirgo) achieve their design sensi- 
tivities it will transform GW astronomy from a search for 
the first detection, into a tool to explore many different 
astrophysical and cosmological phenomena. 

Gravitational waves directly encode the parameters of 
the emitting system, including the luminosity distance 
Dl and the redshifted masses. Simultaneous measure- 
ments of the redshift and the luminosity distance would 
allow gravitational waves to be used as standard can- 
dles, probing the cosmic distance ladder and allowing 
for measurements of cosmological parameters [3, fl5| . 
However, the redshift and the intrinsic masses for point- 
mass objects can not be individually determined from 
gravitational-wave observations alone. Therefore, pre- 
vious attempts to use gravitational waves as standard 
sirens have generally relied on the existence of electro- 
magnetic counterparts which can be used to unambigu- 
ously measure the redshift and break the degeneracy 
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we demonstrate that such counterparts are not necessary 
if the intrinsic mass distribution is sufficiently narrow, as 
may be the case for double NS (DNS) binaries, although 
one can do even better by combining the two approaches. 

We show that it is possible to use the statistics from a 
catalog of gravitational-wave observations of inspiraling 
DNS systems to simultaneously determine the underlying 
cosmological parameters and the NS mass distribution. 
A given cosmological model determines the redshift as 
a function of luminosity distance, making it possible to 
extract the intrinsic mass of a system from a measure- 
ment oi and the redshifted mass. This permits us 
to statistically constrain the Hubble constant and the 
NS mass distribution via a Bayesian formalism, using 
only GW data. A narrower intrinsic NS mass distribu- 
tion will more effectively penalize any model parameters 
which arc offset from the true values. We investigate how 
the precision with which we can recover the underlying 
parameters scales with the number of detections and the 
values of the intrinsic parameters themselves. 

For the majority of our analysis, we do not consider 
difficult-to-detect electromagnetic (EM) counterparts to 
the GW detections, which were relied on in previous anal- 
yses, e.g., [l7j | . Nor do we consider tidal coupling cor- 



rections to the phase evolution of DNS inspiral signals, 
which break the degeneracy between mass parameters 
and redshift to probe the distance- redshift relation [l9j . 
but which only enter at the fifth post-Newtonian order 
and will likely be very difficult to measure with Advanced 
LIGO. Rather, we rely on measurements of the redshifted 
chirp mass, which is expected to be the best-determined 
parameter, and the luminosity distance. This approach 
was introduced by Markovic in |20j, where the author 
extracted candidate source redshifts from the redshifted 
chirp mass using a constant intrinsic chirp mass (this 
is later extended to include some spread around the as- 
sumed intrinsic value). Chernoff and Finn explored this 
technique in [2l[ , which was elaborated upon by Finn in 
, where he suggested using the distribution of signal- 
to-noise ratios and chirp masses to probe cosmological 
parameters. In this paper, we use up-to-date cosmol- 
ogy, mass-distribution models, expectations for detector 
sensitivity and parameter measurement accuracies to in- 
vestigate the precision with which the Hubble constant, 
and NS mass distribution parameters, could be measured 
by the advanced GW detector network. 

The paper is organized as follows. In Sec.|Hl we present 
a simplified analytical calculation and derive scaling laws. 
Section |HT] describes the assumptions made in creating a 
catalog of sources, including a discussion of the DNS sys- 
tem properties we can deduce from a gravitational wave 
signal, as well as neutron-star mass distributions and 
merger rates. Section HVl details the theoretical machin- 
ery for analyzing a catalog of detected DNS systems and 
the details of our implementation of this analysis. We 
describe our results in Sec. [Vj in which we illustrate the 
possibility of probing the Hubble constant and neutron- 
star mass distribution via GW data, and conclude in Sec. 
IVII with a summary and discussion of future work. 



II. ANALYTICAL MODEL 

Here, we present a simplified analytical model that we 
use to show the feasibility of our idea and to derive the 
main scaling relationships. We provide additional justi- 
fication for the various assumptions made in this model 
later in the paper. 

The network of advanced detectors will be sensitive 
to gravitational waves from NS-NS binaries only at rel- 
atively low redshifts, z < 0.15 (see Sec. IIII A[) . At such 
low redshifts, the Hubble law is nearly linear, so that to 
lowest order, we can write the Hubble constant as (see 
Section HIEI 
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Therefore, we expect that the uncertainty in the extrap- 
olation of Hq from redshift and distance measurements 
will scale as 
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The detected neutron-star binaries will yield a catalog 
of sources with measured parameters. These parame- 
ters will include estimates of the redshifted chirp mass 
M. z = (1 + z)M. and luminosity distance Dl- The red- 
shifted chirp mass will be measured very accurately, so we 
can ignore measurement errors for this parameter. How- 
ever, our ability to extract the redshift of an individual 
source from the redshifted chirp mass will depend on the 
narrowness of the intrinsic chirp mass distribution, 
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where the last approximation follows from the fact that 
z <C 1. On the other hand, the luminosity distance is 
estimated directly from the gravitational- wave signal, but 
with a significant error that is inversely proportional to 
the signal-to-noise-ratio (SNR) of the detection. 

Existing binary pulsar measurements suggest that the 
chirp-mass distribution may be fairly narrow, um 
0.06Afo (see Sec. MI Dj) . Meanwhile, for the most dis- 
tant sources at the threshold of detectability, z s» 0.15 
and \SD L \/D L » 0.3 (see SecHE). Therefore, the first 
term in Eq. ([2J is generally larger than the second term 
(though they become comparable for the most distant 
sources), and the intrinsic spread in the chirp mass dom- 
inates as the source of error. 

The errors described above were for a single detec- 
tion, but, as usual, both sources of uncertainty are re- 
duced with more detections as 1/y/N, where N is the 
total number of detected binaries. In principle, we could 
worry whether a few very precise measurements domi- 
nate over the remaining ~iV, affecting the overall 1/y/N 
scaling. The term (a M /M)/z is larger than \SD L \/D L , 
so the best measurements will be those where the for- 
mer term is minimized. The spread in the intrinsic chirp 
mass aj^/Ai is independent of the SNR. Thus, we will 
learn the most from measurements at high z, even though 
these will have a worse uncertainty in Dl (the SNR 
scales inversely with Dl). Therefore, somewhat counter- 
intuitively, the low SNR observations will be most infor- 
mative. However, since the detections are roughly dis- 
tributed uniformly in the volume in which the detector 
is sensitive, we expect half of all detections to be within 
^20% of the most distant detection: therefore, we do ex- 



pect a oc 1/y/N scaling in SHo/ Hq. 

Using the values quoted above, for N ~ 100 detections, 
we may expect that it will be possible to extract the 
Hubble constant with an uncertainty of ~5%. We carry 



1 This scaling holds whenever the number of detections is in- 
creased, either because the merger rate is higher or because data 
are taken for longer. On the other hand, if the number of detec- 
tions increases because the detectors become more sensitive, the 
distance or redshift to the furthest detection will also increase, 
scaling with iV 1 / 3 . In that case, as long as the first term in 
Eq. J2} is still dominant, the overall improvement in SHq/ Ho 
scales as 1/7V 5 / 6 . 



out a rigorous analysis below, and find that the results 
of our simplistic model are accurate to within a factor of 
-2 (see Sec. El). 



III. SOURCE CATALOG 



A. System properties from the gravitational 
waveform 



In the following analysis we consider an advanced 
global network detecting the gravitational radiation from 
an inspiraling double-neutron-star system. The wave- 
form for such an inspiral has a distinctive signature. Such 
systems are denoted chirping binaries due to the charac- 
teristic frequency increase up to the point of coalescence. 

We use the formalism of [22J to describe the response 
of an interferometric detector to the gravitational radia- 
tion from an inspiraling binary. The detector response is 
a function of the system's redshifted mass, the luminos- 
ity distance to the source and the relative orientation of 
the source and detector. This relative orientation is de- 
scribed by four angles. Two of them (6, <f>) are spherical 
polar coordinates describing the angular position of the 
source relative to the detector. The remaining two (t, ip) 
describe the orientation of the binary orbit with respect 
to the observer's line of sight [22j |. 

In the quadrupolar approximation, the dependence 
of the detector response, h(t), on these four angles is 
completely encapsulated in one variable, 0, through the 
equation 



-6(7r/) 2 / 3 cos[ X + $(t)], fort<T, 

for t > T, 



(4) 



where / is the GW frequency, x 1S a constant phase, 
<!>(£) is the signal's phase, and T is taken as the time of 
coalescence. M z = (l+z)A4 is the redshifted chirp mass, 
while M. encodes an accurately measurable combination 
of the neutron star masses, 
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M 



mim2 



(mi + m,2) 



{mi + m 2 ). 



(5) 



Analysis of the gravitational- wave phase evolution yields 
errors on the deduced redshifted chirp mass which vary 
according to the waveform family being used. Regardless, 
the precision is expected to be extremely high, with a 
characteristic error of ^0.04% 2 1231] . 



9 is defined as 



6 = 2[F 2 (1 + cos 2 if 



4^ cos 2 t] 1 ' 2 , 



(6) 



2 M 2 can be determined from the strain signal in one interferom- 
eter through the phase evolution. 
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where < O < 4, and 

F + = — (1 + cos 2 9) cos 2(f) cos 2ip — cos 9 sin 20 sin 2ip, 

F x = - (1 + cos 2 9) cos 2(f) sin 2ip + cos sin 2<f> cos 2t/>. 

(7) 

The luminosity distance Dl is encoded directly in the 
gravitational waveform; however a single interferometer 
cannot deduce this. The degeneracy in the detector re- 
sponse between O and Dl must be broken and this re- 
quires a network of three or more separated interferom- 
eters to triangulate the sky location Jl4]. A network of 
separated detectors will be sensitive to the different GW 
polarization states, the amplitudes of which have differ- 
ent dependences on the binary orbital inclination angle, 
l. Thus, the degree of elliptical polarization measured 
by a network can constrain i [9j. The interferometers 
comprising a network will be misaligned such that their 
varying responses to an incoming GW can constrain the 
polarization angle, ip. 

Once is constrained, Dl can then be deduced from 
the detector response, giving a typical measurement error 
of ~(300/p)%, where p is the signal-to- noise ratio of the 
detection (e.g., |23l - t25| ). The accuracy with which the 
distance can be measured will depend on the exact net- 
work configuration (for example, the introduction of an 
Australian detection will partially break the inclination- 
distance degeneracy [26j]), but we will use the above as a 
representative value. 

We don't include the impact of detector amplitude cal- 
ibration errors, which could lead to systematic biases in 
distance estimates. Unlike statistical measurement er- 
rors, these biases would not be ameliorated by increasing 
the number of detections. For example, calibration errors 
of order 10%, as estimated for the LIGO S5 search [27| . 
would translate directly into 10% systematic biases in 
Hq estimates. Thus, systematic calibration errors could 
become the limiting factor on the accuracy of measuring 
Ho if they exceed the statistical errors estimated in this 
paper. 

B. Detector characteristics 

For the purposes of creating a catalog of sources for our 
study, we are only interested in determining which bina- 
ries are detectable, and how accurately the parameters 
of these binaries can be estimated. We use the crite- 
rion that the network signal-to-noise ratio, p nc t, must be 
greater than 8 for detection. Actual searches use signifi- 
cantly more complicated detection statistics that depend 
on the network configuration, data quality, and search 
techniques, which might make our assumed detectability 
threshold optimistic. Here, we are interested only in a 
sensible approximation of the detectability criterion. 

The network configuration for the advanced detector 
era is uncertain at present. Possibilities include two 



LIGO 4-km detectors at Hanford and one at Livingston 
(HHL), probably sharing data with a 3- km European 
Virgo detector (HHLV). Alternatively, one of the Hanford 
detectors may be moved to Australia (AHL or AHLV), 
improving the network's parameter-estimation accuracy 
[26|,[28j], while the Japanese detector LCGT and/or an 
Indian detector IndlGO may join the network at a later 
date. 

In the HHL configuration all of the sites are located 
in the United States, such that we may use the ap- 
proximation of assuming the AdLIGO interferometers 
can be used in triple coincidence to constitute a super- 
interferometer. This assumption is motivated by the ori- 
entation of the interferometer arms being approximately 
parallel 1291 . and also has precedents in the literature 
[e.g.. [22T |30||. However, source localization and Dl de- 
termination is very poor in HHL, and would be greatly 
improved by the inclusion of data from Virgo or an Aus- 
tralian detector. 

The single-interferometer approximation is less obvi- 
ously valid for networks with distant, nonaligned detec- 
tors, such as AHL(V) or HHLV. In [31(, the authors com- 
ment that the proposed LIGO- Australia site was chosen 
to be nearly antipodal to the LIGO sites such that all 
three interferometers in the AHL configuration would 
have similar antenna patterns. Furthermore, since the 
same hardware configuration would be used for LIGO- 
Australia and AdLIGO, the noise spectra are expected 
to be similar [32|]. Meanwhile, Virgo does not have the 
same antenna pattern as the LIGO detectors, and the 
Advanced Virgo noise spectrum 33] will be somewhat 
different from the Advanced LIGO spectrum. 

In any case, precise comparisons of the sensitivity of 
different networks depend on assumptions about search 
strategies (e.g., coincident vs fully coherent searches) and 
source distributions (see, e.g., [26|, HJ [34| ) . We therefore 
penalize our super-interferometer assumption in two dif- 
ferent ways. Firstly, we set the network SNR threshold 
to correspond to the expected SNR from three identi- 
cal interferometers, as described below, rather than the 
four interferometers comprising the AHLV or HHLV net- 
works. We further penalize the HHLV network relative to 
the network including the more optimally located LIGO- 
Australia by raising the SNR threshold from 8 to ~10. 
These increases in SNR thresholds have the effect of re- 
stricting the network's reach in luminosity distance or 
redshift; however, similar numbers of detections can be 
achieved by longer observation times. 

With the aforementioned caveats, we proceed with our 
assumption that a global network can be approximated as 
a single super-interferometer. This is to provide a proof 
of principle for the ability of such a network to probe the 
background cosmology and aspects of the source distri- 
bution. We do not anchor our analysis to precise knowl- 
edge of the individual interferometer site locations and 
orientations, but will attempt to correct for any possible 
bias. 

Following [22| (and correcting for a missing square 



5 



root), we can write the matched filtering signal-to- noise 
ratio in a single detector as 
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S/t(/) denotes the detector's noise power spectral den- 
sity (the Fourier transform of the noise auto-correlation 
function) , and 2/ max is the wave frequency at which the 
inspiral detection template ends (35[. The SNR of a de- 
tected system will vary between the individual network 
sites, as a result of the different Sh(fYs and angular de- 
pendencies. The network SNR of a detected system is 
given by the quadrature summation of the individual in- 
terferometer SNRs, 
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We approximate the sensitivity of the super- 
interferometer by assuming 3 identical interferometers 
in the network with the sensitivity of Advanced LIGO, 
such that To. n ot ~ v3ro- Different target noise curves 
for AdLIGO produce different values for rn, which vary 
between ~ 80 — 120 Mpc [32| . We adopt the median 
value of 100 Mpc for a single interferometer, yielding 
?"o,net ~ 176 Mpc for the network. 

The SNR also depends on C(/max), which increases 
monotonically as a function of / max . This factor describes 
the overlap of the signal power with the detector band- 
width [22], which will depend on the wave frequency at 
which the post-Newtonian approximation breaks down, 
and the inspiral ends. It is usual to assume that the in- 
spiral phase terminates when the evolution reaches the 
innermost stable circular orbit (ISCO), whereupon the 
neutron stars merge in less than an orbital period. This 
gives 



f GW 
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l + Z 



1570 Hz /2.8M, 
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where M is the total mass of the binary system [12j . 
/isco also depends directly on the mass ratio p/M (p 
is the system's reduced mass); however this mass asym- 
metry term has a negligible effect on / max for the mass 
range of neutron stars considered here [35L [3(| . 

The maximum binary system mass could conceivably 



be ~4.2M Q . 3 The AdLIGO horizon distance for 1.4M Q - 
1.4M inspirals is ^445 Mpc, which corresponds to 
z ~ 0.1 in the ACDM cosmology. Given that we are 
evaluating different cosmological parameters, we adopt 
z ~ 1 as a generous upper redshift limit to a second- 
generation network's reach. This redshift exceeds the 
reach of AdLIGO in all considered cosmologies 4 and chirp 
masses. With these extreme choices for the variables, 
the orbital frequency at the ISCO, / max , could be as 
low as ^262 Hz. For the latest zero-detuning-high-power 
AdLIGO noise curve [H, C(/max = 262Hz) > 0.98. 
Thus, we feel justified in adopting C(/max) — 1 for the 
ensuing analysis. 

Thus matched filtering, with an SNR threshold of 8, a 
characteristic distance reach of ~176 Mpc and C(/max) — 
1, provides a criterion to determine the detectability of a 
source by our network. 



C. Orientation function, Q 



The angular dependence of the SNR is encapsulated 
within the variable O, which varies between and 4, and 
has a functional form given by Eq. ([5]). From our cat- 
alog of coincident DNS inspiral detections we will use 
only M. z and Dl for each system. The sky location and 
binary orientation can be deduced from the network anal- 
ysis, however we will not explicitly consider them here. 
Without specific vales for the angles (8, 0, t, ip) we can 
still write down the probability density function for O 
[22| . Taking cos 6, <f>/w, cost and ip/n to be uncorre- 
cted and distributed uniformly over the range [—1,1], 
the cumulative probability distribution for O was calcu- 
lated numerically in [37j . The pro bability distribution 
can be accurately approximated 22] by, 



0, 



e) 3 



if < < 4, 
otherwise. 



(12) 



3 Both neutron stars in the binary system would need to have 
masses 2cr above the distribution mean at the maximum consid- 
ered fi and cr, where fi-^s £ [1.0, 1.5]Mq, (Tns £ [0, 0.3]M@. 

4 H £ [0.0,200.0] km s^Mpc" 1 ; U kfi = 0; fi m ,o e [0.0,0.5] 

5 There will be some bias in this approximation, since we are 
assuming each interferometer records the same SNR for each 
event. The fact that the different interferometers are not colo- 
cated means that this may overestimate the number of coincident 
detections. We carry out the analysis here aware of, but choosing 
to ignore, this bias, and in Sec. IV Cl consider raising the network 
SNR threshold, which has the same effect as reducing the char- 
acteristic distance reach of the network. 
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We can use Eq. (jl2|) to evaluate the cumulative distribu- 
tion of 9, 



C e (a 



T e (Q)dG 



1, 

(l+x)(4-x) 4 
256 

0, 



if x < 
if < x < 4 
if x > 4. 

(13) 



D. Mass distribution 

In recent years, the number of cataloged pulsar binary 
systems has increased to the level that the underlying 
neutron-star mass distribution can be probed. There is 
now a concordance across the literature that the neu- 
tron star mass distribution is multimodal, which reflects 
the different evolutionary paths of pulsar binary systems 
38, 39]. However, we are only concerned with neu- 
tron stars in NS-NS systems for this analysis, and their 
distribution appears to be quite narrow. In particular 
39] found that the neutron stars in DNS systems pop- 
ulate a lower mass peak at m ~ 1.37M Q ± 0.042Af Q . 
Meanwhile, [38| restricted their sample of neutron stars 
to those with secure mass measurements and predicted 
that the posterior density for the DNS systems peaked 
at /i NS ~ 1.34M©, cr NS ~ 0.06M Q . 

Population synthesis studies of binary evolution pre- 
dict similarly narrow mass distributions for neutron stars 
in NS-NS binaries (see, e.g., [l(| H3, E3| and references 
therein). Some models predict that the mass of neutron 
stars at formation is bimodal, with peaks around 1.3 and 
1.8 solar masses, and any post-formation mass transfer 
in DNS systems is not expected to change that distribu- 
tion significantly, but the 1.8M mode is anticipated to 
be very rare for DNS systems, with the vast majority of 
merging neutron stars belonging to the 1.3M Q peak 41 1. 
Thus, population synthesis results support the anticipa- 
tion that NS binaries may have a narrow range of masses 
that could be modeled by a Gaussian distribution. 

To lowest order, the GW signal depends on the two 
neutron star masses through the chirp mass, M.. We 
assume that the distribution of individual neutron-star 
masses is normal, as suggested above. For ons <C MnSj 
this should yield an approximately normal distribution 
for the chirp mass as well. 

We carried out ^O(10 5 ) iterations, drawing two ran- 
dom variates from a normal distribution (representing 
the individual neutron-star masses) , and then computing 
M. . We varied the mean and width of the underlying dis- 
tribution within the allowed ranges (Sec. IIIIBI) . Binning 
the M. values, the resulting A4 distribution was found to 
be normal, as expected. 

We postulate a simple ansatz for the relationship be- 
tween the chirp mass distribution parameters and the 
underlying neutron star mass distribution. If X\ and Xi 
are two independent random variates drawn from normal 



distributions 

X X r, 

aXi 



bX 2 ~ N(afii 



X 2 



b cr 2 ). 
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Since the neutron-star mass distribution is symmet- 
ric around the mean (and all neutron star masses are 
~O(1M ) with the values spread over a relatively nar- 
row range), then we can assume a characteristic value 
for the pre-factor in Eq. is the value taken when both 
masses are equal i.e. ^(0.25) 3 / 5 . The chirp mass distri- 
bution should then be approximately normal 

M~N(» c ,a 2 c ), 

with mean and standard deviation 

fi c « 2(0.25) 3 / 5 /iNS, a c « V2(0.25) 3/5 a NS - (15) 

where /ins and 0ns are the mean and standard deviation 
of the underlying neutron-star mass distribution, respec- 
tively. 

The accuracy of such an ansatz depends upon the size 
of mass asymmetries which could arise in a DNS binary 
system. We investigated the percentage offset between 
the actual distribution parameters (deduced from least- 
squares fitting to the sample number-density distribu- 
tion) and the ansatz parameters, for a few values of /^ns 
and <tns- The largest offset of the ansatz parameters 
from the true chirp mass distribution was on the or- 
der of a few percent (~2.5% for fi c , and ^3.5% for a c 
when //NS = 1-OMq, ctns = 0.3M©), and the agreement 
improved with a narrower underlying neutron star mass 
distribution. For the case of ctns ~ 0.05M©, the agree- 
ment was ~0.1% for fi c and < 0.1% for a c . In the case of 
<7ns ~ 0.15Mq, the agreement was within a percent for 
both parameters. The sign of these offsets indicates that 

,,true - ,, model —true _modcl 

M c < M c and a c > a c 

Given that the literature indicates an underlying 
neutron-star mass distribution in DNS systems with 
cns ^ 0.15M Q , we anticipate that Eq. (fTS")) will be ap- 
propriate for generating data sets and we use this in the 
ensuing analysis. The assumption throughout is that for 
the volume of the Universe probed by our global network, 
the neutron-star mass distribution does not change. 

The observed data will tell us how wide the intrinsic 
chirp mass distribution is in reality. If it is wider than an- 
ticipated, we may not be able to measure Hq as precisely 
as we find here, but we will know this from the observa- 
tions. In principle, we could still be systematically biased 
if the mass distribution turned out to be significantly 
non-Gaussian, since we are assuming a Gaussian model. 
However, it will be fairly obvious if the mass distribution 
is significantly non-Gaussian (e.g., has a non- negligible 
secondary peak around 1.8 Mq), since redshift could only 
introduce a ~10% spread in the very precise redshiftcd 
chirp mass measurements for detectors that are sensitive 
to z ~ 0.1. In such a case, we would not attempt to fit 
the data to a Gaussian model for the intrinsic chirp mass 
distribution. 
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TABLE I: A compilation of NS-NS merger rate densities in 
various forms from Tables II, III and IV in [l2| . The first column 
gives the units. The second, third and fourth columns denote the 
plausible pessimistic, likely, and plausible optimistic merger rates 
extrapolated from the observed sample of Galactic binary neutron 
stars [44H . The fifth column denotes the upp er rate limit deduced 
from the rate of Type Ib/Ic supernovae [45l| . 



Source 




Rrc 


^high 


Rmax 


NS-NS (MWEG~ 1 Myr _1 ) 


1 


100 


1000 


4000 


NS-NS (L^Myr" 1 ) 


0.6 


60 


600 


2000 


NS-NS (Mpc^Myr" 1 ) 


0.01 


1 


10 


50 



E. DNS binary merger rate density, h(z) 

We assume that merging DNS systems are distributed 
homogeneously and isotropically. The total number of 
these systems that will be detected by the global network 
depends on the intrinsic rate of coalescing binary systems 
per comoving volume. We require some knowledge of this 
in order to generate our mock data sets. Any sort of 
redshift evolution of this quantity (as a result of star- 
formation rate evolution etc.) can be factorised out [22j . 
such that 



fi(t) = 



d 2 N 
dt„dV r 



= h(z) = n £(z), 



(16) 



where N is the number of coalescing systems, t c is proper 
time, V c is comoving volume, and ho represents the local 
merger-rate density. 

We will consider an evolving merger-rate density, such 
that, 



1 



az 



l + 2z, 



for z < 1, 



42 



(17) 
to the 



which is motivated by a piecewise linear fit 
merger-rate evolution deduced from the UV-luminosity- 
inferred star- formation-rate history (43|. 

The appropriate value for fiQ is discussed in detail in 
12]. In that paper, the authors review the range of 
values quoted in the literature for compact binary co- 
alescence rates, i.e. not only NS-NS mergers but NS-BH 
and BH-BH. The binary coalescence rates are quoted 
per Milky Way Equivalent Galaxy (MWEG) and per L w 
(10 10 times the Solar blue-light luminosity, Lb,®), as well 
as per unit comoving volume. In each case, the rates are 
characterized by four values, a "low," "realistic," "high" 
and "maximum" rate, which cover the full range of pub- 
lished estimates. 

The values for the NS-NS merger rate given by [l2[ 
are listed in Table HI The second row of Table U is de- 
rived assuming that coalescence rates are proportional 
to the star-formation rate in nearby spiral galaxies. This 
star-formation rate is crudely estimated from their blue- 
luminosity, and the merger-rate density is deduced via 
the conversion factor of 1.7 L 10 /MWEG 46]. The data 
in the third row is obtained using the conversion factor 
of 0.0198 Lio/Mpc 3 0. 



To convert from merger-rate densities to detection 
rates, [ijj take the product of the merger-rate density 
with the volume of a sphere with radius equal to the vol- 
ume averaged horizon distance. The horizon distance is 
the distance at which an optimally oriented, optimally lo- 
cated binary system of inspiraling 1.4M neutron stars is 
detected with the threshold SNR. This is then averaged 
over all sky locations and binary orientations. 



47T 

N D = n a x — 



D 



horizon 

Mpc 



(2.26)" 



(18) 



where the (1/2.26) factor represents the average over all 
sky locations and binary orientations. 

This gives ~40 detection events per year in AdLIGO 
(using i? ro ), assuming that -Dhorizon = 445 Mpc and all 
neutron stars have a mass of 1.4M . 



F. Cosmological model assumptions 

We assume a flat cosmology, Qk,o = 0, throughout, for 
which the luminosity distance as a function of the radial 
comoving distance is given by 



D L (z) = (l + z)D c {z) = (l + z)D H 



where Djj = c/H (the "Hubble length scale") and 
E(z) = ^Jn m , (l + z) 3 + Q A ,o- 



(19) 



(20) 



In such a cosmology, the redshift derivative of the comov- 
ing volume is given by 



dVc = 4ttD c (zYD h 
dz E(z) 



(21) 



At low redshifts, we can use an approximate simplified 
form for the relationship between redshift and luminos- 
ity distance. Using a Taylor expansion of the comoving 
distance around z — up to 0(z 2 ), and taking the ap- 
propriate positive root, we find D c (z) = Dl(z)/(1 + z) is 
given by 



D c {z) w D c {z = 0) + z 



dD c 



z 2 8 2 D r 



2! dz 2 



z=0 



-.QmfiZ' 



(22) 



Hence, 



D L ^D H 



z + ( 1 - -fi m ,o ) ^ 



Therefore, 



2 (l — -Sl m o) 



/ i | 4 (i - f r> m , ) d l i 



(23) 
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This approximation is very accurate for the range of 
parameters investigated (H Q <G [0,200] km s _1 Mpc _1 , 
^m.o € [0,0.5]), and for Dl < 1 Gpc (which is comfort- 
ably beyond the reach of AdLIGO for NS-NS binaries). 
In this parameter range, the largest offset of this ap- 
proximation from a full redshift root-finding algorithm is 
~4.6%, at a luminosity distance of 1 Gpc. 



G. Distribution of detectable DNS systems 

The two system properties we will use in our analysis 
are the redshifted chirp mass, A4 Z , and the luminosity 
distance, Dl- Only systems with an SNR greater than 
threshold will be detected. Thus, we must include SNR 
selection effects in the calculation for the number of de- 
tections. We can write down the distribution of the num- 
ber of events per year with M , z and 6 [H, |3(| , 



d A N 



dV c h(z) 
dtdOdzdM dz (1 + z) 



V{M)V B {@), (24) 



where t is the time measured in the observer's frame, such 
that the 1/(1 + z) factor accounts for the redshifting of 
the merger rate [3Cj|. 

We convert this to a distribution in M. z , D L and p 
using, 



d^N 



dtdpdD L dM z 



dM 


dM 


dM 


dM- 


dD L 


dp 


dz 


dz 


dz 


dM- 


dD L 


dp 


ae 


d& 


96 


dM z 


dD L 


dp 



d N 



dtdQdzdM 



(25) 



We use the definitions of the variables in Sec. IIII Al and 
IIIIBI to evaluate the Jacobian matrix determinant. The 
redshift is only a function of Dl (in a given cosmology) ; 
the intrinsic chirp mass, M, is the redshifted chirp mass 
divided by (1 + z) (again the redshift is a function of 
Dl); 6 is a function of M z , Dl and p according to Eq. 

®. The (1,3) component (dM/dp = (dM/dp)\ Mi D ) 

is zero because we are differentiating intrinsic chirp mass 
(a function of redshifted chirp mass and distance) with 
respect to SNR, but keeping distance and redshifted chirp 
mass constant. If these variables are held constant then 
the derivative must be zero. Similar considerations of 
which variables are held constant in the partial deriva- 
tives are used to evaluate the remaining elements of the 
matrix. Hence, 



(l+z) 



5 e 

6 M z 



We note that, 



M dz 
(l+z) dD L 
dz 
dD L 

D L 



1 



dz e 



(1 + z) 8D L p ' 



(26) 



v p (p)Sp = r e (e)se, 



which gives, 



V p (p\M z ,D L )=Ve(e)^ 



d P Mz,D 
p£l (1.2M q 
8 r " 



P 



5/6' 



M z 



D, 



1.2M, 



5/6 



8r V M z 
such that we finally obtain, 

d 4 N _ 1 dz dV c h(z) 



(27) 



dtdpdD L dM z (1 + z) ODl dz (1 + z) 

xV(M\z)x -Pe(e)- 
P 

" v ' 

V p (p\Mz,D L ) 



4irD c {z) 2 D H 



h(z) 



D c (z)E(z) + D H (l + z) (l + z) 2 
M 2 



x V 



l + z 



D L )V P (p\M z ,D L ). 



(28) 



We may not necessarily care about the specific SNR 
of a detection; rather only that a system with M. z 
and Dl has SNR above threshold (and is thus de- 
tectable). Fortunately the SNR only enters Eq. (|28|) 
through V p {p\M. z ,Dl), such that we can simply inte- 
grate over this term and apply Eq. (|13p . 



POO POO 

/ V p (p\M z ,D L )dp= V e (Q)de = C e (x), 

J On J x 



Po D L (l.2M & 
where, x = — 



ro 



M z 



5/6 

(29) 



In this case, Eq. (|28[) is modified to give, 
d 3 N 



dtdD L dM z 



P>P0 



4ttD c (z) 2 D h 



h(z) 



D c {z)E(z) + D H (l + z) (l + z) 2 



x V 



M 2 



l + z 



d l a 



Po£l (\J2Mq 
8 r V M z 



5/6' 



(30) 



To calculate the number of detected systems (given a set 
of cosmological and NS mass distribution parameters, 7^) 
we integrate over this distribution, which is equivalent to 
integrating over the distribution of events with redshift 

and chirp mass, i.e. N p — T x f™ J™ ( rft £^ ) dzdM, 

where T is the duration of the observation run. 
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TABLE II: A summary of the WMAP 7-year observations. 
The data from Column 1 is from Table 3 of 48], containing 
parameters derived from fitting models to WMAP data only. 
Column 2 contains the derived parameters from Table 8 of [49(1 , 
where the values result from a six-parameter flat ACDM model fit 
to WMAP+BAO+SNe data. 



Parameter 



WMAP only WMAP+BAO+SNe 



H I (km s^Mpc -1 ) 

&A,0 



71.0 + 2.5 
0.0449 ± 0.0028 
0.222 + 0.026 
0.734 ± 0.029 



70.4±i;| 
0.0456 ± 0.0016 
0.227 ± 0.014 
0.728t°-™ 



year, which compares well to the 40 events found in p3 



IV. ANALYSIS METHODOLOGY 

We will use Bayesian analysis techniques to simulta- 
neously compute posterior distribution functions on the 
mean and standard deviation of the intrinsic NS mass 
distribution (in DNS systems) and the cosmological pa- 
rameters given a catalog of simulated sources with mea- 
sured redshifted chirp masses and luminosity distances. 



H. Creating mock catalogs of DNS binary 
inspiraling systems 

The model parameter space we investigate is the 5D 
space of [Ho ,o, a] with a flat cosmology as- 

sumed. To generate a catalog of events, we choose a set 
of reference parameters, motivated by previous analysis 
in the literature. The seven-year WMAP observations 
gave the cosmological parameters in Table HU For our ref- 
erence cosmology, we adopt H = 70.4 km s _1 Mpc -1 , 
^m.o = 0.27 and f^A.o = 0.73. The parameters of the 
neutron-star mass distribution were discussed earlier, but 
as reference we use uns = 1.35M and ctns = 0.06Af©. 
The merger-rate density was also discussed earlier, and 
we take a = 2.0 and no = 10 -6 Mpc _3 yr _1 as the ref- 
erence. Later, we will investigate how the results change 
if the width of the NS mass distribution is as large as 
O.13M0, as indicated by the predictive density estimate 

of m. 

These reference parameters are used to calculate an 
expected number of events, 6 and the number of observed 
events is drawn from a Poisson distribution (assuming 
each binary system is independent of all others) with 
that mean. Monte-Carlo acceptance/rejection sampling 
is used to draw random redshifts and chirp masses from 
the distribution in Eq. (|24|) for each of the N a events. 
The Dl and M. z are then computed from the sampled 
M. and z. 

With a reference rate of no = 10~ 6 Mpc _3 yr _1 and 
a constant merger-rate density, we estimate that there 
should be ~90 yr -1 detections, whilst taking into account 
merger-rate evolution using Eq. (fl"7f boosts this to ^100 
yr -1 . These numbers are for a network SNR threshold of 
8. If we ignore merger-rate evolution and raise the SNR 
threshold to 10 (to represent an AdVirgo-HHL network 
for which the coincident detection rate is roughly halved 
relative to the HHL-only network) we get ~45 events in 1 



The observation time, T, is assumed to be 1 year (but the ex- 
pected number of detections simply scales linearly with time) 
and a network acting as a super-interferometer with ro,not — 176 
Mpc is also assumed. 



A. Bayesian analysis using Markov Chain Monte 
Carlo techniques 

Bayes' theorem states that the inferred posterior prob- 
ability distribution of the parameters ~jl based on a hy- 
pothesis model %, and given data D is given by 



C(D\~jt,H)Tr(~jt\H) 
P{D\H) 



(31) 



where C{D\~jl ,1-L) is the likelihood (the probability of 
measuring the data, given a model with parameters ft), 
(jt\H) is the prior (any constraints already existing 
on the model parameters) and finally p(D\H) is the ev- 
idence (this is important in model selection, but in the 
subsequent analysis in this paper can be ignored as a 
normalization constant). 

In this analysis, the data in Eq. pip is not from a 
single source, but rather from a set of sources, and we 
want to use it to constrain certain aspects of the source 
distribution, as well as the background cosmology. The 
uncertainty arises from the fact that any model cannot 
predict the exact events we will see, but rather an as- 
trophysical rate of events that gives rise to the observed 
events. The probability distribution for the set of events 
will be discussed in Sec. IIVBI 

To compute the posterior on the model parameters, 
we use Markov Chain Monte Carlo (MCMC) techniques 
since they provide an efficient way to explore the full pa- 
rameter space. An initial point, xq, is drawn from the 
prior distribution and then at each subsequent iteration, 
i, a new point, l/, is drawn from a proposal distribution, 
q(y | x) (uniform in all cases, covering the range of pa- 
rameter investigation). The Metropolis-Hastings ratio is 
then evaluated, 



R = 



n(l})£(D\l},H)q&\T) 



(32) 



A random sample is drawn from a uniform distribution, 
u G U[0, 1], and if u < R the move to the new point is 
accepted, so that we set Xi+{ =~y- If u > R, the move 
is rejected and we set Xi + \ = ~x%. If R > 1 the move is 
always accepted, however if R < 1 there is still a chance 
that the move will be accepted. 
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The MCMC samples can be used to carry out integrals 
over the posterior, e.g. 



J fC#)p&\D, H)d^ = 1 /(*$)■ (33) 



The ID marginalized posterior probability distributions 
in individual model parameters can be obtained by bin- 
ning the chain samples in that parameter. 



B. Modelling the likelihood 

1. Expressing the likelihood 

We use a theoretical framework similar to that of 50] . 
The data are assumed to be a catalog of events for which 
redshifted chirp mass, A4 Z , and luminosity distance, Dl 
have been estimated. These two parameters for the 
events can be used to probe the underlying cosmology 
and neutron-star mass distribution. In this analysis, we 
focus on what we can learn about the Hubble constant, 
Hq, the Gaussian mean of the (DNS system) neutron- 
star mass distribution, /ins j and the Gaussian half- width, 
ctns- We could also include the present-day matter den- 
sity, Q m 0) however we expect that this will not be well 
constrained due to the low luminosity distances of the 
sources. We could also include the gradient parameter, 
a, describing the redshift evolution of the merger-rate 
density. 

The measurement errors were discussed earlier (Sec. 
IIII A|) and we will account for these later. For the first 
analysis we assume that the observable properties of in- 
dividual binaries are measured exactly. 

We consider first a binned analysis. We divide the 
parameter space of M. z and Dl into bins, such that the 
data is the number of events measured in a particular 
range of redshifted chirp mass and luminosity distance. 
Each binary system can be modeled as independent of 
all other systems, so that within a given galaxy we can 
model the number of inspirals that occur within a certain 
time as a Poisson process, with DNS binaries merging at 
a particular rate (e.g., dHHH). The mean of the Poisson 
process will be equal to the model-dependent rate times 
the observation time, and the actual number of inspirals 
occurring in the galaxy is a random- variate drawn from 
the Poisson distribution. 

A bin in the space of system properties may contain 
events from several galaxies, but these galaxies will be- 
have independently and the number of recorded detec- 
tions in a given bin will then be a Poisson process, with 
a mean equal to the model-dependent expected number 
of detections in that bin [Hfj. The data can be written 
as a vector of numbers in labelled bins in the 2D space of 
system properties, i.e. ~n = (711,712, • ■ -,nx), where X is 
the number of bins. Therefore, the likelihood of record- 
ing data D under model H (with model parameters ~ft ) is 



the product of the individual Poisson probabilities for de- 
tecting rii events in a bin, i, where the expected (model- 
dependent) number of detections is /•;("/£). For the 
bin, 



: ill 



p(n % \~]l ,U) 



n,,\ 



(34) 



and so the likelihood of the cataloged detections is, 

^i7^)=n (r ' (7?)) T' ,7?1 - <*> 



In this work, we take the continuum limit of Eq. (|35p . 
In this case, the number of events in each infinitesimal 
bin is either or 1. Every infinitesimal bin contributes a 
factor of e~ Ti ^\ whilst the remaining terms in Eq. (|35[) 
evaluate to 1 for empty bins, and ri(jl) for full bins. The 



product of the exponential factors gives e **, where 
is the number of DNS inspiral detections predicted by the 
model, with parameters 7* • The continuum likelihood of 
a catalog of discrete events is therefore 



-N, 



N 

]>(^) 



(36) 



where A = {A~i, A2, . . ., A~v*} is the vector of measured 

system properties, with A, = (A4 Z , Dl)% for system i, and 

N Q is the number of detected systems. Finally, r(Aj I'jt) is 
the rate of events with properties M. z and Dl, evaluated 
for the i th detection under model parameters ~jt , which 
is given by Eq. pi)) . 



2. Marginalizing over ho 

We may also modify the likelihood calculation to 
marginalize over the poorly constrained merger-rate den- 
sity, no- 7 This quantity is so poorly known (see Table HJ, 
that it is preferable to use a new statistic that does not 
rely on the local merger-rate density, by integrating the 
likelihood given in Eq. (|36p over this quantity, 



— j roo —t 

C(A\~tt,H)= £(A\-j},H)dn . (37) 
Jo 



The expected number of detections described in Sec. IIII Gl 

can be expressed as, 



N„ 



n x 



1 dM z dD L , 



7 A similar technique was used in [53j . where the total number of 
events predicted by the model is marginalized over. 
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where, 



4ttD c (z) 2 D h l + az 

' D c (z)E(z) + D H (l + z) (1 + z) 2 



x V 



M 2 



l + z 



D L )C e 



and, 



r$\~j}) = h x J i: 



Po^l /1.2M G 



r V ^« 



5/6' 



(38) 



(39) 



where li is the integrand evaluated for the i th system's 
properties. Thus, 

C(A\-jt,H) 




dh 1 x 



i=l 



The integral, J / I dMzdD^, depends on the underlying 
model parameters, ~jt , through I, but it does not depend 
on riQ. Therefore, defining 

7 = no X J J I dM z dDL = n x 5. 

We note that no € [0, oo], hence 7 <E [0, 00]. Therefore, 



d'y 



%\ ■ < 

SJ S 



N 



7^ e -^ 7 ) xJ-^+DJJl,. (41) 



i=l 



independent of 7^ 



We will verify in the analysis that this new likelihood pro- 
duces results completely consistent with the case where 
exact knowledge of the merger-rate density is assumed. 
We note that we did not include a prior on no in the 
above, which is equivalent to using a flat prior for no £ 
[0, oo]. This reflects our current lack of knowledge of the 
intrinsic merger rate, although such a prior is not nor- 
malizable. We could implement a normalizable prior by 
adding a cut-off, but this cut-off should be set sufficiently 
high that it will not influence the posterior and therefore 
the result will be equivalent to the above. 



C. Calculating the posterior probability 

The likelihood statistic C is used to marginalize over 
the poorly constrained local merger-rate density. We 



use a weakly informative prior on the model parame- 
ters, so that it doesn't prejudice our analysis. As a prior 
on /xns we take a normal distribution with parameters 
fi = 1.35M Q , a = 0.13Mq. This is motivated by the 
posterior predictive density estimate for a neutron star 
in a DNS binary system given in [38j . We take a prior on 
a that is a normal distribution, centred at 2.0 with a a of 
0.5. Uniform priors were used for the other parameters. 

We made sure that the size of the sampled param- 
eter space was large enough to fully sample the pos- 
terior distribution, so that we could investigate how 
well gravitational-wave observations alone could con- 
strain the cosmology and neutron-star mass distribution. 
The parameter ranges were H £ [0, 200] km s _1 Mpc _1 , 
O m ,0 £ [0,0.5], /ins € [1.0,1.5]M Q , a NS £ [O,O.3]M and 
a £ [0.0,5.0]. 

To calculate C for a given point in the model pa- 
rameter space we must compute the number of detec- 
tions predicted by those model parameters (N^), and 
we need to calculate z(Dl) for that model so that A4 
can be evaluated. For the sake of computational effi- 
ciency, some approximations are used. We have verified 
that our results are insensitive to these approximations. 
Our approximation for z(Dl) was described in Sec. lIIIFl 
We also used an analytic ansatz to calculate the model- 
dependent expected number of detections, based on fac- 
torizing the contributions from different model param- 
eters. The agreement between this ansatz and the full 
integrated model number is excellent, with the biggest 
discrepancy being < 3% of the true value. This allows 
a direct calculation of without a multi-dimensional 
integration for each point in parameter space. 



V. RESULTS & ANALYSIS 

For our first analysis, we will assume that A4 Z and Dl 
for each individual merger are measured perfectly by our 
observations, so that the A4 Z and Dl recorded for the 
events are the true values. This represents the best case 
of what we could learn from GW observations. Later, 
we will consider how the accuracy of the reconstructed 
model parameters is affected by including measurement 
errors on the recorded event properties. 



A. Posterior recovery 

We carried out the analysis discussed in Sec. IIVI for 
a data set calculated from the reference model given in 
Sec. IIII HI We found that ^10 6 — 10 7 samples were nec- 
essary for the MCMC analysis to recover the underlying 
posterior distributions. 

We found that O mi o and a were not constrained by the 
observations, but their inclusion in the parameter space 
did not affect our ability to recover the other parameters. 
For this reason, we kept them in the analysis, but all 
remaining results will be marginalized over these model 
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FIG. 2: Recovered 2D posterior distribution in Ho and /^ns 
space, showing a correlation between these recovered parameters. 
The model parameter values used to generate the data are the 
reference values. There appears to be negligible correlation 
between cxns with Hq or /xns- 



parameters. Given the low redshift range that a second- 
generation network is sensitive to, it is not surprising 
that the matter-density and merger-rate evolution were 
not constrained. 

The recovered ID posterior distributions in the other 
parameters are shown in Fig. [T] for a typical realization 
of the set of observed events. We have verified that these 
marginalized distributions are consistent with those ob- 
tained when exact knowledge of the intrinsic no is as- 
sumed. We found that the ID posterior distributions for 
.Ho, ln(<7Ns) and /zns were well fit by Gaussian distribu- 
tions of the form Aexp {—{x — /i) 2 /2a 2 ). These best-fit 
Gaussians are also shown in the Figure. Although the 
distributions do not peak at the model parameters used 
to generate the data, those values are consistent with the 
mean and width of the recovered distributions. 

In Fig. [2] we show the corresponding 2D posterior dis- 
tribution in Ho and /iNS parameter space. We see that 
a correlation exists between these parameters. Given a 
cataloged Dl value, a low value of Hq will imply a low 



model-dependent redshift. When this redshift is used 
to compute M from M z , we calculate a large value of 
the chirp mass, which implies a chirp mass distribution 
(and hence a neutron-star mass distribution) centered at 
larger values, ons simply encodes the width of the mass 
distribution around the mean, so on average it should 
have no effect on Hq and /ins calculations and indeed 
we found that ctns showed no correlation with the other 
model parameters. 

It is clear from Fig.[T]that the parameters of the Gaus- 
sian fits provide a useful way to characterize the recov- 
ered distributions. We can then describe the recovered 
distributions in terms of two best-fit parameters i.e. the 
Gaussian mean, /i, and Gaussian half-width, a. 

B. Random spread of best-fit parameters 

1. No errors in data catalog 

To explore the spread in the best-fit parameters of the 
recovered posteriors over different realizations of the data 
catalog, we generated 100 different realizations, keeping 
the intrinsic parameter values the same for each. 

In each case, we fit a Gaussian to the ID posteriors 
and record the mean, /x, standard deviation, er, and the 
"error" in the mean. This last quantity is the number 
of standard deviations that the mean is offset from the 
intrinsic value, i.e. AX = (/x — X)/a, where X is the 
value of the parameter used to generate the catalog [5(| ■ 
A ±2ct offset encloses ~95% of the Gaussian probability 
distribution, so we would reasonably expect most of the 
realizations to lie within this range. 

Figure [3] shows the distributions of the Gaussian-fit 
standard deviations and "errors" for Hq, ln(c7Ns) and 
/iNS over 100 realizations of the AdLIGO-network data 
catalog. The distribution of the Gaussian-fit means for 
each parameter roughly resemble their respective poste- 
riors, and the distribution of Gaussian standard devia- 
tions also appears approximately Gaussian. As we would 
have hoped, most of the realizations have a best-fit mean 
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FIG. 3: Distribution of the Gaussian-fit standard deviations (top) and "errors" (bottom) of the recovered posteriors over 100 
realizations, for Hg (left), ln(<7Ns) (center) and ^tNS (right). More details are given in the text. 



which is offset from the intrinsic value by less than 2a. 
As with the Gaussian-fit parameters, the "error" distri- 
bution is approximately Gaussian and centered around 
i.e. centered around the intrinsic value. 

The most useful quantity here is the standard deviation 
of the reconstructed posterior distribution, as it charac- 
terizes how well we will be able to constrain the model pa- 
rameters. The distribution over 100 realizations displays 
the typical range of this "measurement accuracy." Thus, 
ignoring measurement errors in the data, and with refer- 
ence parameters used to generate the catalog, we could 
conceivably determine Hq, ctns and p^s to an accuracy 
of -±10 km s^Mpc -1 , -±0.004M Q , 8 and ~±O.O12M 
respectively. 



2. Including & accounting for errors 



As discussed in Sec. MI Al the system properties of each 
event in the catalog will include some error arising from 
instrumental noise. The data for each event will actually 
be in the form of posterior probability density functions 
(PDFs) for the properties, where previously we have as- 
sumed these are ^-functions at the true values. We repeat 



the analysis assuming uncertainty in the source proper- 
ties. We can include errors in the system properties in 
the data generation stage, by choosing the recorded val- 
ues from a Gaussian distribution centered on the true 
value, with a standard deviation of 0.04% for A4 Z and 
(300/p)% for D L , where p is the SNR of the detected 
event. 



When we included errors in the data generation, but 
did not account for them in the analysis, we found that 
the model parameter posterior distributions were on av- 
erage biased toward lower values of Hq, with biases also 
present in the p^s and 0ns distributions. When the er- 
rors are added, systems will move both to lower and to 
higher values of the luminosity distance. However, as 
we discussed in Sec. II, the sources at greatest distance 
have the most influence on our ability to measure the cos- 
mology. We would therefore expect the sources shifted 
to greater distances to have most impact on the cosmo- 
logical parameter estimation, biasing us toward smaller 
values of Hq, as we found. 



Evaluated using <5(ctns) = ""NS x < 5( m (< T Ns))i taking ctns to be 
the reference value and a typical error in In (ons) of 0.072. 



However, we can account for these errors in the analy- 
sis, by modifying the previous likelihood in Eq. Q36p (5G | 
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FIG. 4: A comparison of the best-fit a distributions over fOO realizations, between the case of no errors present in the data catalog, 
and the case of errors applied to system properties in the catalog. We have attempted to compensate for the errors in the data, (i), (iii) 
and (v) show the best-fit a distributions when no errors are applied to the system properties in the data catalog, (ii), (iv) and (vi) show 
the best-fit a distributions when the received data has errors. 



to 



-N„ 



N 



p I n = s 



d k >^d k )^---d k xZ, (42) 



where, in our case, each system is associated with two 



cataloged properties such that k = 2, and s is the de- 
tector output, which is a combination of N a signals, hi, 
and noise, n. This is as an integral over all possible 
values of the source parameters that are consistent with 
the data. The first term inside the square bracket is the 
computed posterior PDF for the detected population of 
sources. Typical LIGO/ Virgo DNS inspiral detections 
last only a few seconds, whilst AdLIGO/ AdVirgo inspi- 
ral detections may be in-band for several minutes. Re- 
gardless, these detections should be uncorrelated, with 
independent parameter estimates [54] . and so this first 
term reduces to the product of the posterior PDFs for 
each detection. 

If the posterior PDF for a given source has been 
obtained via MCMC techniques, then the integral in 
Eq. (|42l) may be computed by summing over the chain 
samples. Thus, errors may be accounted for by making 
the following replacement in Eq. (|36l) . 



-Kj) 



(43) 



where Mi is the number of points in the chain for the i th 
source's PDF, and X^ is the j th element of the discrete 
chain representing this PDF. This technique does not 
assume a specific form for the PDF, and can be used 
in the case of multimodal distributions. 9 



In this analysis, we include errors on Dl only, as those 
on the redshifted chirp mass M. z are very small and can 
be ignored. (The uncertainty in the redshift estimate, 
which dominates the uncertainty in Hq as discussed in 
Sec. |nl arises from the width of the intrinsic chirp-mass 
distribution.) We represent the Dl posterior PDF for 
each source by a chain of 75 points, drawn from a normal 
distribution with standard deviation a — (3/p)£>L, and 
a mean equal to the value in the data catalog, which in 
this analysis, as discussed earlier, includes an error to 
offset it from the true value. Whilst we adopt a simple 
Gaussian Dl posterior PDF, the methodology we use 
here to account for errors is not reliant on the specific 
form of the PDF. 

Using this analysis, we found that the bias in the pos- 
terior means for Hq, ctns and Mns was corrected. In Fig. 
5] we show a comparison of the best-fit a distributions 
for each of the parameters when measurement errors are 
included (and accounted for), compared to the case in 
which they are ignored. It is clear that the presence of 
measurement errors decreases the measurement precision 
that we can achieve. However, the distributions overlap 
in all cases, and the peak of the error distributions is 
shifted only ^20% higher. 

These errors only cause a shift in the measurement 
precision, so that we can ignore errors in the cataloged 
properties, with the knowledge that a full analysis would 
produce broadly the same results, but with ^20% worse 
precision. The presence of errors (when accounted for) 
should therefore not affect our general conclusions about 
what a second-generation global network will be able to 
tell us about the underlying cosmological and source pa- 



9 Multimodal distributions may result from partial degeneracies 



with other waveform parameters [54| . such as the angular vari- 
ables encapsulated in 0. Examples of this are shown in |17| . 
where the sky position of a detected system is pinned down, and 
the degeneracy between the inclination angle, t, and can lead 
to multimodal posteriors for Dl which skew the peak to higher 
distances than the intrinsic value. 
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rameters. 



C. Dependence on number of observed events 

The next question we will explore is how the measure- 
ment accuracy of the parameters depends on the number 
of cataloged events. This can be answered by chang- 
ing the local merger-rate density, no, or the observation 
time, T, whilst keeping the other model parameters fixed. 
We analyzed catalogs with different values of ho around 
the previously used realistic value (2.5xl0 -7 , 5.0xl0~ 7 , 
l.OxlCT 6 and 2.0xl0~ 6 Mpc _3 yr- 1 ), using 10 realiza- 
tions in each case. 

In Fig. [3] we show the standard deviation of the 
recorded posterior distribution versus the number of cat- 
aloged events for each realization of each fiQ. The distri- 
butions are well fit by a function of the form, 



1 



a x oc 



(44) 



which one might expect; we have a population of N a 
events which we are using to statistically constrain a pa- 
rameter, so we expect that the root-mean-squared error 
on the parameter should scale as ^l/V^Vo- The points 
and solid lines are the data and best-fit curves when we 
ignore measurement errors in the data generation, whilst 
the dashed lines are best-fit curves to data where we ac- 
count for measurement errors, as in the previous section. 

Table IIIII shows the percentage fractional accuracy to 
which we could measure each model parameter, in both 
the case that we ignore errors and when we account for 
them. The range of local merger-rate densities reflects 
the quoted values in [l2| , and the means of the posterior 
distributions are taken as the reference values. 

The number of detected events will also depend on 
the SNR threshold, po.net In practice, the network 
thresholds required for detection are often higher than ~8 
because a network performs more trials of the same data 
and is sensitive to both gravitational wave polarizations 
simultaneously. The result of increasing the threshold 
SNR to 10 is to approximately halve the detection rate. If 
the expected detection rate is ^100 yr _1 in the po, nc t = 8 
case, this becomes ~50 yr _1 in the po, nc t = 10 case. 

This halving of the detection rate is expected since, 



Vc,eff(po 



10) 



Vc,eS (P0,net — 8) 



10 



= 0.512, 



(45) 



where V Cje s is the effective comoving volume to which the 
network is sensitive. We can achieve the same number 
of detections at higher SNR thresholds by increasing the 
observation time. Using a higher network SNR threshold 
is equivalent to assuming a lower characteristic distance 
reach for the network. By increasing po to 10, we cut the 
detection rate in half, which is roughly the decrease in the 
number of coincident detections when we shift from the 



HHL to HHLV network [34(. A network SNR threshold 
of 12 reduces the detection rate to ~30 yr^ 1 . 

To investigate the dependence of the Hq measurement 
accuracy on the characteristic distance reach of the net- 
work (a prediction of our scaling arguments), we com- 
puted 10 realizations at each of 10 different network SNR 
thresholds, ranging from 6 to 15. The detection rates 
were kept the same at each SNR threshold by rescaling 
h Q . The reference values were po,nct = 8 and r 0inc t = 176 
Mpc, as used previously. At each po, n ctj a weighted mean 
of the Gaussian-fit half-widths of the parameter posteri- 
ors was calculated, with error bars determined by the 
maximum and minimum half-widths out of the 10 real- 
izations. The results for Hq are shown in Figure [75] The 
fit favors a (l/Vo,net) relationship, as expected from scal- 
ing arguments. There appears to be no effect on the 
measurement accuracy of the NS mass distribution pa- 
rameters. No measurement errors were included on ei- 
ther the recorded Dl or Ai z values, and the detection 
rate was fixed in this analysis, so it is unsurprising that 
the measurement precision of the NS mass distribution 
parameters is unaffected by the reach of the network. In 
this particular investigation, given that the total num- 
ber of events is unchanged, and therefore the number 
of masses to which we fit the NS mass distribution is 
unchanged, we do not expect the precision of the fit to 
change either. 

This indicates that the measurement accuracies of cns 
and /xns quoted in this paper will be achievable at differ- 
ent po,net and n^net by scaling the observation times, or 
if the Universe has a different ho than expected. How- 
ever, the measurement accuracy of Hq is also linked to 
the characteristic distance reach of the network. 



D. Dependence of measurement accuracy on 
intrinsic parameters 

It is also interesting to investigate how the constraints 
on the parameters of the underlying distributions depend 
on the parameters used to construct the distribution. 
This was done by generating 10 data realizations at each 
of 24 different combinations of the intrinsic pns an d 0ns ■ 
The intrinsic values of Hq, f2 TOi o and a were fixed at their 
reference values. The recorded measurement precision for 
a given intrinsic parameter combination was the weighted 
mean of this value over 10 realizations. Figure [7] shows 
the results of this analysis. One can see that the mea- 
surement precision depends on the width of the intrinsic 
NS mass distribution. An increase in the intrinsic ctns 
by a factor of 6 leads to a reduction in the measurement 
accuracy on Hq and pns by a factor of ^6, but only leads 
to a modest 10% reduction of the measurement accuracy 
for ln(er NS ). 

The improvement of the measurement accuracy with 
a narrower intrinsic DNS mass distribution is a key re- 
sult. In order to constrain the Hubble constant to within 
^±10% with ^100 observations, we require the Gaussian 
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FIG. 5: The "measurement accuracy" of each parameter (represented by the standard deviation of the Gaussian fit to the posterior) 
plotted against the number of observed events, N . The intrinsic parameters are kept fixed whilst the local merger-rate density, no, is 
scaled up and down. The number of observed events scales linearly with the observation time and the local merger-rate density, such 
that the same result is achieved for twice the local merger-rate density if the observation time is halved. We see that a 1/V ' N relation is 
favored. The points and solid lines correspond to the case when we ignore errors, where the curves have gradients 108±2 km s _1 Mpc _1 , 
0.737±0.001 and O.131±O.OO2Af0 respectively. The dashed lines are best-fit curves for the same analysis, but with measurement errors 
included and accounted for, for which the gradients are 136 km s _1 Mpc _1 , 0.917 and O.152M0, respectively. 



TABLE III: The best-fit curves to the plots in Fig. [5] are used to compute the percentage measurement precision of the model 
parameters. The local merger rates match the range quoted in [T3 |. where in our analysis, no = 1.0 Mpc - 3 Myr — 1 gives ~100 detections 
in 1 year (at a network SNR threshold of 8). In each case, the mean of the posterior distribution is taken at the reference value. 
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FIG. 6: The variation of measurement accuracy with 
instrumental distance reach is shown. Each point represents the 
weighted mean of the Ho measurement accuracy from 10 
realizations at a particular ro net, where the error bars show the 
maximum and minimum values of cr out of the 10 values. All 
other parameters are at their reference values, and the total 
number of detections is scaled to 100. The reference distance 
reach is ro, n ct ~ 176 Mpc. The curve is a (l/r"o, n et) fit to the 
data, with gradient 10.8 ± 0.2 km s^Mpc -1 . 



half-width of the DNS mass distribution to be smaller 
than 0.04Mq. The explanation for this is that we esti- 
mate the system chirp mass, M, by dividing the red- 
shifted chirp mass, A4 Z , by (1 + z), where the z is 
model-dependent (having been calculated from Dl with 
given cosmological parameters). Thus, a narrower NS 
mass distribution will more effectively penalize model pa- 
rameters which deviate from the intrinsic values. For 
ons = O.13M [38j . an accuracy of ^±10% on iJ would 
require ~O(1000) detections. 

The dependence of the measurement precision on /xns 
is not very clear from the left and right panels, but the 
effect on <tns is evident in Fig. Ojfb)] Varying the intrinsic 
//ns from 1.33M Q to 1.39M Q provides a ^5 — 10% gain in 
In (cns) precision. The variation of the expected number 
of detections with ctns is less than one event, whilst for 
/ljns it is more significant. So, all the posterior fit a values 
were scaled to the average number of detections for a 
given (Uns- This varies by ~15 detections over the range 
of /*ns investigated. 

To explain the improvement in measurement precision 
with larger values of //ns, we note Eq. ([5]). We see that 
a larger mass distribution mean will, on average, imply 
larger individual NS masses. For a fixed SNR thresh- 
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FIG. 7: The variation of the weighted mean (over 10 realizations) of the Gaussian-fit standard deviations with the parameters of the 
underlying NS mass distribution. All other parameters are fixed at their reference values. The variation of the expected number of 
detections with <r N g is less than one event, whilst for /^ns it is more significant. Thus all the posterior fit a values are scaled to the 
average number of detections for a given fi^s e.g. for /xns = 1.33Mq this average event number is ~95, whilst for 1.39-Mq it is ~110, in 
1 year. 



old, this allows detections to be made from larger Dl 
values, thereby raising the effective comoving volume to 
which the network is sensitive. This raises the number 
of detections, and hence the parameter measurement ac- 
curacy. Rescaling all the measurement accuracies to 100 
events confirms that this is the dominant effect, as the 
different /^ns curves in Fig. [7] then overlap. Factorizing 
out the dependence on N a also confirms that the varia- 
tion of the measurement accuracy with the width of the 
underlying NS mass distribution is a real feature. 

Repeating the above analysis for fixed ^nSj but with 
different combinations of Hq and uns confirms the varia- 
tion of precision with ctns- However, there appears to be 
no strong dependence on H as it is varied by ±10 km 
s~ x Mpc _1 around the reference value. 10 

E. Complementing GW data with GRB redshift 
data 

In the above, we have assumed only the GW observa- 
tions are available. However, if the redshift of the system 
is somehow known, then the background cosmological pa- 
rameters can be directly probed using the luminosity dis- 
tance, Dl, measured from the GWs (l4l |. 

Merging compact-object binary systems, such as NS- 
NS or NS-BH, are leading candidates for the explana- 
tion of short-duration gamma ray bursts (SGRB). SGRB 
events are among the most luminous explosions in the 
universe, releasing in less than one second the energy 
emitted by our Galaxy over one year [55j . and involv- 
ing intense outflows of gamma rays. There is therefore a 
good chance that EM counterparts to GW-detected DNS 
mergers will be observed. There is strong evidence that 
the emission from GRBs is not isotropic [56l458l |. which 



The reference value is well constrained by WMAP+BAO+SNe 
analysis [49ll , 



may be due to the formation of relativistic jets in these 
systems (55|. The redshift can be determined from the 
longer- wavelength SGRB afterglow (59|. 

Therefore, we may only observe the event electromag- 
netically if we happen to lie within the cone of the ra- 
diative outflow, whilst we should be able to detect the 
gravitational wave signal from any DNS merger within 
the AdLIGO-network horizon. We denote the beaming 
fraction, /, for SGRBs with a double-jet of opening angle 
0j by [Mil, 

'-£-—(!)■ < 46) 

If we assume the population of SGRBs is randomly ori- 
ented on the sky, that their progenitors are all NS-NS 
mergers, that we will detect all SGRBs that are beamed 
toward us, and that the required SNR for a GW detection 
is independent of the existence of a counterpart, then the 
beaming fraction, /, is also the fraction of DNS inspiral- 
ing systems for which we would be able to gather red- 
shift data. In practice, GW searches that are triggered 
by electromagnetic observations of SGRBs would have a 
greater distance reach than blind analyses, which would 
tend to increase the fraction of counterparts. However, 
gamma-ray telescopes operating in the advanced detec- 
tor era might not have 100% sky coverage, which would 
tend to reduce the fraction of counterparts. In addition, 
even with an SGRB counterpart we might not be able to 
determine the redshift, as this requires observation of an 
afterglow. However, all of the GW sources for advanced 
detectors will be at low redshift, for which the chances 
of measuring the redshift are significantly higher. In the 
following, when we refer to the beaming fraction, we will 
mean the fraction of GW detections with electromagnet- 
ically determined redshifts, which will be similar to the 
intrinsic beaming fraction, but not exactly the same for 
the reasons just described. 

We performed a simple analysis to see how the mea- 
surement accuracy would improve if some fraction, /, of 
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FIG. 8: Plots of the measurement precision versus the beaming fraction of SGRBs, /. Approximately fN events in the catalog were 
denoted as SGRBs, while the remaining ~(1 — f)N events in the catalog were assumed to be GW-only events. A single data catalog 
was used repeatedly, with larger and larger fractions of it assumed to be observable as SGRBs. The data was generated with the 
reference parameters. The fitted curves are of the form aexp(— by/x) + c, where (a, b, c) = (9.04, 3.64, 1.59) and (0.00932, 4.46, 0.00439) 
for the Hq and /^ns precisions respectively. The corresponding plot for ctims shows no obvious trend. 



redshift data was available. A single data set was gener- 
ated with the redshift, luminosity distance and redshifted 
chirp mass of each event recorded (with reference intrin- 
sic parameters). The measured M. z and Dz, were drawn 
from Gaussian distributions centered at the true value, 
as described previously. However, as before the small er- 
ror on M z was ignored. When only a gravitational- wave 
signal is available the system properties are analyzed as 
previously, with the measured values assumed to be the 
true values. 



(1411) splits into two components, 



If an event is included in the SGRB fraction (with an 
associated redshift), then the likelihood is the product 
of the GW likelihood with the redshift posterior PDF, 
which we take to be a delta function (since spectroscopic 
redshift determination will be much more precise than 
GW determinations of Dl). The Dl posterior was taken 
to be a Gaussian, centered around the measured value, 
and with a standard deviation of 30% of this distance. 
This percentage error is a worst case, corresponding to 
a detection near the threshold SNR, i.e. (300//9)% for 
p = 10. Using a constant percentage of the distance as 
the width of the Dl posterior is pessimistic, since closer 
events will be measured with greater accuracy. Integrat- 
ing over redshift in Eq. (|4"2j) picks out the value of the 
integrand at the true system redshift, since the redshift 
posterior is a delta function. Thus, the product in Eq. 
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The identification of the SGRB as an EM counterpart 
will vastly improve sky localization of the source, helping 
to beat down the degeneracies in the GW observations 
between D^ and the inclination angle, t. A fuller in- 
vestigation could also consider a prior on the inclination 
angle, given that the source is an SGRB with a colli- 
mated outflow [It], HH, E3] , and that emission has been 
observed. This would further help to improve the mea- 
surement precision of the luminosity distance. 

If / = 1, we find the precision for Hq is ~2.0 km 
s _1 Mpc _1 , compared to ~11.0 km s _1 Mpc _1 when / = 
0. The results are shown in Fig. HI along with fits to the 
data of the form a exp(-by^x) + c. The important result 
here is that the accuracy with which we are able to con- 
strain Hq and (i^s improves markedly with the beaming 
fraction. This is to be expected, since by recording z 
and M z we know exactly what the intrinsic chirp mass, 
A4, of the system is. The high accuracy of the redshift 
measurements restricts the space of model parameters 
through the Gaussian factor in Eq. (|4"7) . The same plot 
for ctns shows no trend at all. This may be because the 
measurement accuracy of ctns is most strongly linked to 
the number of cataloged events, rather than whether we 
include extra system information. 
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This analysis could be sensitive to the errors we include 
in the data catalog, since the normal distribution in the 
left square-bracket of Eq. (|47| will favor model parame- 
ters, ~jt, such that — DL{zj,~f}). Dj) is the mean 
of the Dl posterior PDF for the j th event, which may be 
skewed away from the true value. However, the intrinsic 
values were always consistent with the mean and width 
of the recovered posteriors, so this does not seem to be a 
significant problem. 

The SGRB jet opening angle is poorly constrained by 
observations. In (6l| . the authors quote the inverse beam- 
ing fraction to be in the range, 1 <C / _1 < 100, giving 
/ > 10~ 2 or a jet opening angle dj > 16°, which is consis- 
tent with theoretical constraints on the jet half-opening 
angle 55] . Such models permit the jet half-opening an- 
gle to be as lar ge a s 30°, for which the beaming fraction 
becomes ~0.13 [55)]. This would allow Hq and ^ns to be 
measured with a precision more than twice that of their 
GW-only values (see Figure [3]). 

In p7| . the authors performed an analysis on multiple 
DNS inspirals detections in the AdLIGO- Virgo network 
with associated EM signatures. They assumed the sky 
location of the sources were known, and that Dl and z 
were measured, so that they directly probed the distance- 
redshift relation. With 4 SGRBs they predicted Hq could 
be measured with a fractional error of ~13%, improving 
to —5% for 15 events. With / = 1, and scaling the mea- 
surement precision as \/^/Wl, we find 4 SGRBs gives 
~12.5% precision, whilst 15 gives ~6.5% precision. Fig- 
ure [8] indicates an Hq precision of ~5% when / = 0.15; 
thus the second square bracket on the right of Eq. (j47|) 
slightly improves the measurement accuracy of Hq com- 
pared to the first square bracket alone. These results are 
dependent on the modelled Dl errors, but are broadly 
consistent with 17]. In contrast, we expect we can con- 
strain Hq to within ^±15% using ^100 GW events, with 
no EM signatures recorded for any of the GW detections. 
15 SGRB events out of the ~100 GW events requires 
a beaming fraction of ~0.15, which is rather optimistic 
given the current constraints on the jet opening angle. 
However this could conceivably be achieved over obser- 
vation times longer than one year; additionally, the de- 
tection of an electromagnetic transient could allow the 
sensitivity volume to be increased in a triggered search. 

In this section, we have not considered the possibility 
of redshift determination of the DNS inspiraling system 
via its association with a host galaxy. This could prove 
difficult in practice, since the sky error box is sufficiently 
large as to contain many candidate galaxies. In [62| the 
authors comment that over 100 galaxies can be found in a 
typical LIGO /Virgo GW signal error box at a distance of 
100 Mpc. However, in the same work they introduced a 
ranking statistic which successfully imaged the true host 
of a simulated GW signal ^93% of the time, if 5 wide- field 
images were taken. The caveat here is that this statistic 
has only been tested out to 100 Mpc, since comprehen- 
sive galaxy catalogs are lacking beyond this range. The 
catalog completeness is not 100% at 100 Mpc, and even if 



more distant, complete catalogs were available, the num- 
ber of potential host galaxies in an AdLIGO/ AdVirgo sky 
error box would be much greater. Dl determination via 
network analysis may help to restrict the redshift range 
of these searches, but this is an area in need of future 
attention. 

A novel method was proposed in [l8[ in the context 
of LISA EMRI detections. In that case, instead of pre- 
cisely identifying the host galaxy of a GW detection (and 
thus the redshift of the source) , the value of Hq was av- 
eraged over all galaxies present in LISA's sky error box. 
Each galaxy in the box was weighted equally, and the 
chosen host galaxy was not included in the likelihood 
calculation to take into account the fact that the true 
host galaxy may not even be visible in available catalogs. 
They showed that sub-percent accuracies on Hq would 
be possible if 20 or more EMRI events are detected to 
z < 0.5. This method has recently been investigated in 
the context of DNS inspirals in the advanced detector era, 
where a precision of a few percent on Hq was claimed to 
be possible with 50 detections (63| . 



VI. CONCLUSIONS & FUTURE WORK 

We have explored the capability of an advanced global 
network of GW interferometers, such as the AHL (Aus- 
tralia, Hanford, Livingston) or HHLV (Hanford, Hanford, 
Livingston, Virgo) configurations, to probe aspects of the 
background cosmology and the nature of the neutron-star 
mass distribution (for NSs in DNS systems). Current rate 
estimates suggest these systems could be a strong candi- 
date for the first direct GW detection. With the reach 
of the advanced detectors, it may be possible to pro- 
duce catalogs of tens of these systems along with their 
associated properties over the first few years of advanced 
detector operation. 

We used a Bayesian theoretical framework to assess 
the posterior probability of the cosmological parameters 
and the mean and standard deviation of the NS mass 
distribution. Catalogs of DNS system mergers were gen- 
erated, comprising the system redshifted chirp mass, A4 Z 
and luminosity distance, Dl, from which we endeavoured 
to statistically constrain the underlying parameters. 

We simulated catalogs of 100 detected binaries (cor- 
responding to a few years of observation for a local 
merger-rate density of 10 -6 Mpc~ 3 yr _1 12]) for refer- 
ence parameters Hq — 70.4 km s _1 Mpc _1 , £! m ,o = 0.27, 
^ns = 1.35A/ , ons = 0.06M Q , a — 2.0 (where a is the 
gradient of the redshift evolution of the NS-NS merger- 
rate density) . With such catalogs of detections we found 
it should be possible to measure the Hubble constant, 
as well as the mean and half-width of the DNS Gaussian 
mass distribution. Hq should be constrained to ^±10 km 
s^Mpc" 1 , ln(cr NS ) to -±0.07 and ^ N s to -±0.012M Q . 
As a result of the restricted cosmological reach of second- 
generation detectors, £l m ,o and a cannot be constrained 
by such observations. This is because the different cosmo- 
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logical density parameters do not significantly affect low 
z luminosity distances, and low z sources will not char- 
acterize the redshift evolution of the merger-rate density. 

The measurement accuracy was characterized by the 
width of a Gaussian fit to the recovered posterior distri- 
butions. We also attempted to account for measurement 
errors in the data catalog and found that taking errors 
into account would slightly broaden the recovered pa- 
rameter distributions, but only by ~20%. This can be 
compensated for by longer observation times. 

Keeping the intrinsic parameters fixed, and scaling the 
merger-rate density (or the observation time) allowed us 
to investigate how this precision varied with the number 
of cataloged events. We found that precisions varied as 

— 1/2 

N for all three parameters. We also investigated 
the effect of changing the network SNR threshold, which 
has the same effect as reducing the distance reach of the 
network. Scaling the local merger-rate densities to give 
equal numbers of detections was enough to achieve the 
same precision on the NS mass distribution parameters, 
but the uncertainty in measuring Ho also scales inversely 
with the distance reach of the network. 

We also checked how the values of the intrinsic param- 
eters themselves affected our ability to constrain them. 
Varying Hq over a range of reasonable values had little 
impact on the measurement precision, but the effect of 
ctns wa s considerable. Changing <tns from 0.12M Q to 
O.O2Af led to a factor of ^6 increase in the precision 
on Hq and /zns> but a modest ~10% improvement on 
In (ctns ) • Our key result is that for Ho to be constrained 
to within ~±10% using ^100 events (with the intrinsic 
Hq and mean of the DNS mass distribution fixed at their 
reference values), then the half-width of the intrinsic DNS 
mass distribution would have to be less than 0.04M Q . 

Finally, considering that NS-NS and NS-BH merger 
events are leading candidates for the progenitors of short- 
duration gamma-ray bursts [55l . [60| , we investigated how 
the measurement precision would improve if redshift data 
were available for some fraction of the catalog. The red- 
shift could be deduced from the afterglow of the SGRB 
or from the closest projected galaxy. The fraction of GW 
detections that have observable EM counterparts will de- 
pend on the opening angle of the SGRB jets. The most 
recent GR-MHD simulations permit a half-opening an- 
gle of 30°, for which the maximum fraction of the DNS 
inspiraling systems that could have an observable EM is 
^0.13 [55] (this fraction could be further increased by the 
greater sensitivity of GW searches triggered on EM tran- 
sients). This would permit a significant improvement on 
the measurement precision of Hq and /zns to more than 
double their GW-only precisions. There appears to be 
no effect on the measurement precision of ctns- 

Our results were based on a single-interferometer for- 
malism to describe the global network, assuming that 
LIGO-Australia would be nearly antipodal to the U.S 
sites and have identical sensitivity. There is no differ- 
ence in the number of expected detections between the 
HHL and AHL configurations, although a slight improve- 



ment in the detection efficiency is expected for the HHLV 
network [24] . We can penalize the number of coincident 
detections made by all interferometers by raising the net- 
work SNR threshold from 8 to 10. This cuts the detection 
rate in half, but this can be compensated for by longer 
observation times. 

We have shown the significant potential for a network 
of second-generation detectors to provide an indepen- 
dent measurement of the Hubble constant, and to de- 
termine the neutron-star mass distribution for those NSs 
found in DNS systems. Even more powerful constraints 
should be possible with the Einstein Telescope (ET), a 
proposed third generation ground-based interferometer 
with an arm- length of 10 km [64|. 

ET will be sensitive to sources out to z ~ 2 for DNS 
inspirals, with the expected number of detections in one 
year being ^O(10 5 — 10 6 ) [64|, [65J. When we compare 
this to AdLIGO's 445 Mpc reach, giving ^100 network 
detections, we see the clear improvement ET will offer. 
With such a large reach and detection rate we anticipate a 
much greater measurement precision on Ho, as well as the 
other parameters discussed in this paper. For the model 
parameters at the reference values, and using our analysis 
to extrapolate for a conservative ET detection rate gives 
~0.5%, ^0.3% and ^0.03% precision on iJ , ctns and 
/Uns respectively. Furthermore, the cosmological reach of 
ET may permit Q m ,o and Slk,o to be constrained, with 
consequences even for probing the dark energy equation 
of state parameter, w j6q.l67j. Preliminary results for our 
future ET analysis has constrained ft m .o to ^±30% with 
100 events (where, for this preliminary study, we used 
the same methodology as in the present paper, but with 
the characteristic distance reach modified to account for 
the sensitivity curve of the early ET design study, ET- 
B [68]), and scaling this for a conservative ET detection 
rate of 10 5 yr" 1 gives ~±0.9%. The ability to detect 
z > 2 events may provide an opportunity to measure the 
evolution of the DNS merger-rate density, which will shed 
light on the evolution of the star-formation rate. Some 
of the techniques used in this paper would have to be 
adapted for any ET analysis, e.g., the approximation to 
deduce the redshift from the luminosity distance would 
have to be replaced by the full root-finding algorithm. 
This coupled with the huge number of catalogued events 
would lead to longer computation times. 

In this analysis we have considered a global second- 
generation GW-intcrfcrometer network. The improve- 
ment offered by a Southern Hemisphere gravitational- 
wave detector would be significant for sky localization 
(though only moderate for distance estimates), but this 
may not be realized. However, even with the HHLV net- 
work, we will still be able to place constraints on the 
underlying model parameters by overcoming the reduced 
coincident detection rate with a longer duration network 
science run. For now, we have shown that if a global net- 
work is successful in detecting populations of inspiraling 
DNS systems, then gravitational wave astronomy can be- 
gin to place independent and interesting constraints on 
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Hq, as well as the neutron-star mass distribution. This 
will be a step toward using gravitational- wave astronomy 
for precision astrophysics. 
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